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Abstract A Bayesian approach to calibrating period- 
luminosity (PL) relations has substantial benefits over 
generic least-squares fits. In particular, the Bayesian 
approach takes into account the full prior distribution 
of the model parameters, such as the a priori distances, 
and refits these parameters as part of the process of 
settling on the most highly-constrained final fit. Ad- 
ditionally, the Bayesian approach can naturally ingest 
data from multiple wavebands and simultaneously fit 
the parameters of PL relations for each waveband in 
a procedure that constrains the parameter posterior 
distributions so as to minimize the scatter of the fi- 
nal fits appropriately in all wavebands. Here we de- 
scribe the generalized approach to Bayesian model fit- 
ting and then specialize to a detailed description of ap- 
plying Bayesian linear model fitting to the mid-infrared 
PL relations of RR Lyrae variable stars. For this exam- 
ple application we quantify the improvement afforded 
by using a Bayesian model fit. We also compare dis- 
tances previously predicted in our example application 
to recently published parallax distances measured with 
the Hubble Space Telescope and find their agreement to 
be a vindication of our methodology. Our intent with 
this article is to spread awareness of the benefits and 
applicability of this Bayesian approach and encourage 
future PL relation investigations to consider employing 
this powerful analysis method. 
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1 Introduction 

The period-luminosity (PL) relations of pulsating vari- 
able stars — typically variables of types RR Lyrae 
(RRL), Cepheid, and Mira — are invaluable tools for 
constructing the rung of the distance ladder that con- 
nects the Milky Way to other nearby galaxies, extend- 
ing to '--^ 5 Mpc. Recent applications of this distance 
measurement technique using Cepheids have success- 
fully mated Cepheid distances to SNe la host galax- 
ies and constrained the Hubble Constant, Hq, to 3.3% 
(Riess et al. 2011). The authors have recently derived 
mid-infrared PL relations for RRL variables (Klein 



et al.|[20TT |, and demonstrated their potential to serve 
as important distance indicators for the Large Magel- 
lanic Cloud. Additionally, continuing studies of Mi- 



ras (Whitelock et al. 2008) confirm their potential to 



provide accurate distances even beyond the reach of 
Cepheids. 

The accuracy and precision of any distance measure- 
ment made using the PL relation of a variable star, or 
any population of variable stars within a distant sys- 
tem, is dominated by the uncertainty of the locally cal- 
ibrated PL relation. The general method is to fit a PL 
relation to the variables for which trigonometric par- 



allax measurements are available (Feast and Catchpole 



1997[ ). For more than the past decade only Hippar- 



cos (original catalog published as M.A.C. Ferryman & 



ESA (19971, and improved reduction by van Leeuwen 



( 2007 1 ) could provide these required local distance mea- 



surements to a significantly large sample of local stars 
with the accuracy necessary. More recently, the Hubble 
Space Telescope Fine Guidance Sensor has been used 
to provide higher accuracy parallax measurements for 
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nine Ccphcids ([Benedict et al.|2007[ ) and five RRL vari- 
ables (Benedict et al. 2011 1. In the coining decade, with 



the launch of the ESA's astrometry mission Gaia, the 
sample size of potential PL relation calibrators and the 
accuracy of their parallax distances will be significantly 
augmented ( Cacciari]|2009 1 . 

PL relations are typically calibrated using straight- 
forward, simple, frequentist statistical techniques (see. 



for example, SoUima et al. (2006) fitting RRL vari 



and Glass and Evans 



ables, Matsunaga et al. (2006) fitting type II Cepheids 



Feast and Gatchpole| ( |1997 l fitting classical Cepheids, 
(|2003|) fitting Miras). At the 



basic level, a PL relation is an equation of the form 
M = alogP -f- /?, where M is the absolute magnitude 
(in a given waveband), a is the slope, P is the period 
(in days), and {3 is the zero point magnitude (which 
itself may be a function of metallicity) . This simple 
linear equation can be reliably fit with the method of 
least squares, and the accuracy of the fit can be assessed 
with the standard deviation of the residuals, a metric 
commonly referred to as the scatter. 

A significant limitation of the least squares regres- 
sion method is that it does not make use of the full 
prior probability distribution of the parallax distances. 
Allowing the distances more flexibility to move within 
their prior distributions translates into posterior dis- 
tances that are more consistent with the fitted PL rela- 
tion and therefore more accurate (on average) than the 
prior distance mean. A second limitation is that this 
traditional method is not easily adapted to fitting PL 
relations of the same variables in different wavebands 
simultaneously. Intelligently combining data from mul- 
tiple wavebands has the potential to produce better fi- 
nal PL relation fits. 

A Bayesian approach for fitting the PL relation 
parameters overcomes these traditional limitations. 



Barnes et al. (2003) discusses in substantial depth the 



application of a Bayesian approach to the Cepheid dis- 
tance scale, using physical pulsational models and ra- 
dial velocity data. In the present work we confine our 
examination to the application of Bayesian methods 
in calibrating, purely plienomenologically, the PL rela- 
tions of pulsating variable stars. We use as our example 



the calibration performed by Klein et al. ( 2011 ). In Sec- 



tion [2] we describe the generalized Bayesian modeling 
approach. In Section[3]we work through the application 
of this Bayesian approach to mid-infrared PL relations 
of RRL variables, as first demonstrated by the authors 
in Klein et al. (2011). We perform a traditional, least 



squares fit to the RRL PL relations and compare with 
the fits from our Bayesian approach in Section |4j Fi- 
nally, in Section [5] we draw conclusions and discuss 
future applications. 



2 Technical Explanation of Bayesian PL 
Relation Fitting 

An excellent and thorough description of Bayesian fit- 



ting of linear models is provided in Gelman et al. ( 2003 ) . 



'Barnes et al. (20031 applies Bayesian analysis to the 



Cepheid distance scale, but does not use linear Bayesian 
model fitting for deriving the Cepheid PL relation. Here 
we review the foundation of the Bayesian approach. 

If we assume that our data, denoted by y, follows 
some pattern or rule or model, as in common in nature, 
and denote the model parameter(s) by 9, then we may 
write the probability of the model being true as p(9) and 
the posterior probability that the model is true given 
our observed data y as p{9\y). The probability p{0) 
is the prior distribution on the model, it is what we 
know before making observations. We can also define 
the likelihood as the probability p{y\9) of observing the 
particular data y conditioned on the model 9. Thus, we 
have the unnormalized Bayes' theorem 



p{9\y)^p{9)p{y\9), 



(1) 



the probability of the model being true given the data is 
proportional to the prior probability of the model times 
the likelihoocQ 

To fit observed data ?; to a model with parameters 9 
we simply evaluate Equation [T] throughout a fine grid 
of values for 9 to create the posterior distribution of 
the model. We can examine this posterior distribution 
(through analysis of repeated random draws from the 
distribution) to find the most likely fit parameters, as 
well as uncertainty in these parameters. Furthermore, 
this posterior distribution reveals the most likely true 
values for the data y, which is of course conditional on 
the prior distributions of the data. 



3 Application to Mid-Infrared RRL Variables 



In Klein et al. (2011) we apply Bayesian model fitting 



to a sample of 76 RRL light curves observed with the 
Wide-Field Infrared Survey Explorer (WISE) ( |Wright| 
et al. 2010). Each of the RRL variables was well- 



observed in three WISE bands (Wl at 3.4, W2 at 4.6, 
and W3 at 12 /im) and their prior distance distributions 
are generated by applying the RRL My-[Fe/H] relation 

^The exact (normalized) form of Bayes' theorem is 

p(e)p(y\e) 



p{d\y) ■ 



p{y) 



(2) 



but for our purposes we implement Equation^ by Monte Carlo 
computer simulation and thus the analytical normalization can 
be ignored. 
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given in Chaboyer ( 1999 ) to their Hipparcos mean flux 



V-band magnitudes, correcting for dust extinction. 

Using the nomenclature of Section [2] we define our 
observed data for each RRL as y = (m, P), the appar- 
ent magnitude and period. The unknown fit parame- 
ters are then 6 = (/i, Mq, a, cr), the distance modulus, 
absolute magnitude zero point, PL relation slope, and 
scatter. We can also define f3 — (/i, Mq, a) so that then 
9 = {/3,a). We put an informative normal prior on each 
11 (from the 1^-band distance estimates) and we put a 
flat prior on everything else. 

We then write our statistical model of the PL rela- 
tionship as 



= A^i + + logio(-Pi/-fo) + , 



(3) 



where fii is the distance modulus for ith RRL, Mq.j is 
the absolute magnitude zero point for the jth WISE 
band at P = Pq, where Pq = 0.50118 day is the mean 
period of the sample, and aj is the slope of the PL re- 
lationship in the jth band. We assume that any extinc- 
tion is negligible in these bands. The error terms are 
independent zero-mean Gaussian random deviates with 
variance (aamij)'^, which describe the intrinsic scatter 
in the m^j about the model, where cr is a free param- 
eter which is an unknown scale factor on the known 
measurement errors, (Tm, ;^|^ We fit the model (Equa- 
tion [S]) using a Bayesian procedure, outlined below and 
explicitly described in Section 4 of Klein et al. (2011 ). 



First, we assume a normal (Gaussian) prior distribu- 
tion on each of the distance moduli with mean /io.i and 
standard deviation ct^q •, as described above. For the 
other parameters in our model (Equation|3|, we assume 
a flat, noninformative prior distribution. 

Our likelihood is normal: 

Pirn, P|/3, a) cc p(m|/3, a) = 7V(X/3, aMia.g{al)), (4) 

where N denotes the multivariate normal distribution. 
Note that p{m,P\9) ex. p{m\9) since P is independent 
of m and doesn't depend on any of the parameters 9. 
In the nomenclature of Section [2] we are solving for 



p{9\y) = p{[}, a\m, P) = p(/3|m, P, (j)p{(T\m, P). 



(5) 



The joint posterior distribution, p{j3,a\m, P), can be 
sampled by first drawing from p((T|m, P) and then, con- 
ditional on that draw, selecting from p{j3\m, P, a). The 
posterior distribution for /?, conditional on the value of 
cr, follows the multivariate normal distribution. 



where /3 is the standard maximum likelihood (weighted 
least squares) solution, /3 = (X^I]~-^X*)~^X^S~^m*. 
Unlike the posterior distribution of /3 (given a), the 
posterior distribution of a, p{(j'^\m, P), does not follow 
a simple conjugate distribution. Instead, the distribu- 
tion follows the form 



p{a^\m,P) 



p{(3)p{<j^)L{m\P,p,a) 
pi(3\m,P,a) 



(7) 



where the prior on (3 is proportional to the informative 
prior on /i, the flat prior on a is p(cr^) cx cr~^, and 
the data likelihood L is the product, over all observed 
magnitudes, of the Gaussian likelihood of the data given 
the model (Equation |3| with all parameters specified. 

We draw samples from our joint posterior distribu- 
tion p{(3, a\m, P) using Equations |6] and [t] in conjunc- 
tion. In practice, we compute]^ p(o^^ |to, P) over a fine 
grid of cr values using Equation [7j and then draw a 
sample of cr from this density. For each sampled cr, 
we subsequently draw a /3 from Equation [6] conditional 
on the drawn cr value. We repeat this process 10,000 
times to characterize the joint posterior distribution. 
Using a large sample from this joint posterior distribu- 
tion, we can compute quantities of interest such as the 
maximum a posteriori slopes and zero points of the PL 
relationship of each WISE band, the intrinsic scatter 
of the data around the PL relationship in each band, 
and the spread in the a posteriori distribution of the 
PL parameters. 



4 Comparison with Traditional Fit 

To demonstrate the improvement in the fit from the 
Bayesian approach, which in turn means an improve- 
ment in the predicted distances resulting from the cal- 
ibrated PL relation, we compare it to a traditional 
least squares regression fit. Using the same prior dis- 
tances (technically, the expectation value of the prior 
distance distributions), period measurements, and ob- 
served WISE mean flux magnitudes we perform a least 
squares fit for the slope a and zero point [3 of the PL 
relation. We calculate the scatter (Icr) as the standard 
deviation of the residuals to the fit. The fit parameters 
and scatter are presented in Table [l] In Figure [l] we 
overplot in red the least squares PL relation fits into 



Figure 5 from Klein et al. (2011) 



(6) 



As expected, the actual parameters of the PL re- 
lations, zero point and slope, are statistically consis- 
tent. Comparison of their Icr scatter, however, illus- 
trates that the Bayesian approach produces a set of 



^The average measurement error, am, is 0.013, 0.013, and 0.045 
mag in Wl, W2, and W3, respectively. 



''Assuming that fi = p. Several iterations show that the posterior 
distribution of a is insensitive to the assumed choice of j3. 
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Table 1 Comparison of least squares (subscript LS) and 
Bayesian (subscript B) fits to the WISE RRL PL relations. 
We calculate the scatter (la) as the standard deviation of 
the residuals to each fit. 



Table 2 Recently published RRL HST parallax distances 
(?) compared with the values which Klein et al. (20111 



band 






lo"LS 


qb 


/9b 


lo-B 


Wl 


-0.420 


-1.675 


0.124 


-0.421 


-1.681 


0.007 


W2 


-0.425 


-1.713 


0.124 


-0.423 


-1.715 


0.007 


W3 


-0.503 


-1.763 


0.149 


-0.493 


-1.688 


0.074 



PL relations with nearly eighteen times lower scatter in 
WISE bands Wl and W2 and two times lower scatter in 
W3. The Bayesian fit's significant reduction in scatter 
primarily results from allowing the posterior distances 
to be fit from within the prior distances' distributions. 
Thus, the scatter of the Bayesian fit more closely ap- 
proaches the true intrinsic scatter. 



5 Conclusions 

Tableland Figure [l] clearly demonstrate the Bayesian 
modeling approach's ability to produce a set of more 
tightly constrained PL relations. This benefit of the 
approach primarily arrises from its taking into account 
the full prior distribution on the distances and allow- 
ing the fit to refine these distance distributions. The 
traditional least squares fitting method instead utilizes 
only the expectation value of these prior distances and 
leaves them unchanged during the fitting procedure, 
hence producing a fit with significantly larger scatter. 

It is natural to inquire as to whether the Bayesian 
model is overly constrained. That is, one must be care- 
ful not to allow the posterior distance distributions to 
diverge significantly from the prior distributions. To 
evaluate this property of the Bayesian fit it is reasonable 
to compute and examine the prior, posterior, and pre- 
diction probability densities, as shown for four exam- 



ples in Figure 7 of Klein et al. (2011 1. Finally, we note 



that the the newest HST parallax measurements for the 
four RRL variables V*RRLyr, V*RZCep, V*SUDra, 
and V*UVOct (?) were not yet available for use as 
distance priors in the Bayesian fit performed in ' Klein| 
et al. ( 2011( ), however the distance posteriors produced 
in that work compare quite well with these newly pub- 
lished parallax distances (see Table [2| . We interpret 
these predictions of the parallax distances as a strong 
vindication of our methodology. 

While it is true that the Bayesian modeling method 
described in this article is not strictly applicable to all 
situations (for example, when fitting data which share 
a common distance prior as is standard for studies of 
pulsating variables in the Large Magellanic Cloud), it 
is nevertheless an invaluable tool for calibrating the 



previously predicted through Bayesian linear model fitting 
of the mid-infrared PL relations. All distances are given in 
parsecs. 



Name 


HST d 


PL Fit d 


Bayesian Prior d 


V*RRLyr 


265 ±9 


253 ± 2 


262 ± 15 


V*RZCcp 


394 ± 30 


381 ± 6 


405 ± 23 


V*SUDra 


704 ± 80 


696 ± 7 


696 ± 40 


V*UVOct 


585 ± 34 


536 ± 4 


553 ± 32 



Galactic PL relations of RRL, Cepheid, and Mira vari- 
ables. This approach also has the potential to improve 
the model fits of other phenomena, such as the calibra- 
tion of other distance indicators (SNe la, planetary neb- 
ulae, tip of the red giant branch, etc). Our intent with 
this article is to spread awareness of the benefits and 
applicability of this Bayesian approach and encourage 
future PL relation investigations to consider employing 
this powerful analysis method. 
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Log(P/P„) Log(P/P„) Log(P/P„) 



Fig. 1 Period-luminosity relations for Wl, W2, and W3 (left to right). Results of the Bayesian fitting procedure are 
shown in black, and the results of the least squares fit are overplotted in red. In each figure, the solid lines show the models' 
predictions of the RRL absolute magnitude, as a function of RRL period. The dashed lines show the ±lcr scatter. The top 
panel of each plot shows the residual spread around the best fit model of each fit procedure. The error bars of the Bayesian 
fit for Wl and W2 are smaller than the diamond markers. 
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